

###################################
d <- read.csv("ChaudoinWilf_FPA_data.csv")

	#1049 respondents in the sample
nrow(d3<-subset(d, v10==1 & consent_tx=="Yes" & is.na(support_12)==F & minutes>=8 & minutes<=32)) #n = 1049


##################
#Figure 1 
##################

install.packages("sciplot")
library(sciplot)

	d<-d3
	attach(d)
	# Summarizes the number 1of respondents that received each treatment
	nbyti<- as.matrix(aggregate(support_12 ~ treatment , data=d, FUN=function(x) sum(!is.na(x)) ))
		n1<-nbyti[1,2]
		n2<-nbyti[2,2]
		n3<-nbyti[3,2]
		n4<-nbyti[4,2]
		n5<-nbyti[5,2]
	# Repeats that step, for the number of respondents by treatment who approved
	d2 <- d[ which(d$support_12 =="1"), ]
	abyti<- as.matrix(aggregate(support_12 ~ treatment, data=d2, FUN=function(x) sum(!is.na(x)) ))
		a1<-abyti[1,2]
		a2<-abyti[2,2]
		a3<-abyti[3,2]
		a4<-abyti[4,2]
		a5<-abyti[5,2]
	# Number of draws from the Beta posterior
		n<-5000
	# Draws from the posterior for each treatment
		pi_1 <- as.matrix(rbeta(n, a1+0.5, n1-a1+0.5))
		pi_2 <- as.matrix(rbeta(n, a2+0.5, n2-a2+0.5))
		pi_3 <- as.matrix(rbeta(n, a3+0.5, n3-a3+0.5))
		pi_4 <- as.matrix(rbeta(n, a4+0.5, n4-a4+0.5))
		pi_5 <- as.matrix(rbeta(n, a5+0.5, n5-a5+0.5))
	# Treatment indicators to make the data for the plot cleaner
			#TREATMENT KEY / CODING
				#	1 = OEP
				#	2 = nia / interdependence
				#	3 = NW / network
				#	4 = null
				#	5 = placebo
		t1 <- as.matrix(rep(1,n, nrow=n, ncol=1))
		t2 <- as.matrix(rep(2,n, nrow=n, ncol=1))
		t3 <- as.matrix(rep(3,n, nrow=n, ncol=1))
		t4 <- as.matrix(rep(4,n, nrow=n, ncol=1))
		t5 <- as.matrix(rep(5,n, nrow=n, ncol=1))
	tis <- as.matrix(rbind(t1,t2,t3,t4,t5))
	pis <- as.matrix(rbind(pi_1,pi_2,pi_3,pi_4,pi_5))

	betas <- as.data.frame(cbind(tis,pis))
	betas$treatmentname = factor(betas$V1, labels = c("Direct","Interdep.","Network","Null","Plac."))

	###
	# Figure 1
	###
	# Approval ratings by treatment group
	#pdf("PercSupport_12_ByTreatment_BetaEsts.pdf")
	lineplot.CI(x.factor = betas$treatmentname, response = betas$V2, type="p", ylab="",xlab="",axes=F, font=2,
			data = betas, main = "", #xlab = "Treatment Group", ylab = "% Supporting Regulation (12)",
			ylim=c(0.35,0.65),
			ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)) )
		mtext("Percent Supporting Regulation", side=3,cex=1.15, font=2, line=2.5)
		axis(2, las=2)
		mtext("Treatment Group", side=1,cex=1.15, font=1, line=2.5)
	#dev.off()


####
#FIGURE 2
#	LOVE OF REGULATION 	(1 and 2)
#	PARTY				(3 and 4)
#	IDEOLOGY			(5 and 6)
###

	####
	# Fig 2, graphs 1 and 2 (of 6) - LOVE OF REGULATION 
	###
		d_reglovl <- d[ which(d$reglove1 == "Too much"), ]	#reglov low == currently, too much reg
			nrow(d_reglovl)#235
		d_reglovh <- d[ which(d$reglove1 == "Too little"), ]	#reglov high == currently, too little reg
			nrow(d_reglovh)#427
	###
	# Support 12
	# Treatment pictures for H/L, H/L
	###

	# Summarizes the number of respondents that received each treatment, by h/l
	nbyti_reglovh <- as.matrix(aggregate(support_12 ~ treatment , data=d_reglovh, FUN=function(x) sum(!is.na(x)) ))
	nbyti_reglovl <- as.matrix(aggregate(support_12 ~ treatment , data=d_reglovl, FUN=function(x) sum(!is.na(x)) ))
		n1_reglovh<-nbyti_reglovh[1,2]
		n2_reglovh<-nbyti_reglovh[2,2]
		n3_reglovh<-nbyti_reglovh[3,2]
		n4_reglovh<-nbyti_reglovh[4,2]
		n5_reglovh<-nbyti_reglovh[5,2]

		n1_reglovl<-nbyti_reglovl[1,2]
		n2_reglovl<-nbyti_reglovl[2,2]
		n3_reglovl<-nbyti_reglovl[3,2]
		n4_reglovl<-nbyti_reglovl[4,2]
		n5_reglovl<-nbyti_reglovl[5,2]
	# Repeats that step, for the number of respondents by treatment who approved
	d2_reglovh <- d_reglovh[ which(d_reglovh$support_12 =="1"), ]
	d2_reglovl <- d_reglovl[ which(d_reglovl$support_12 =="1"), ]

	abyti_reglovh<- as.matrix(aggregate(support_12 ~ treatment, data=d2_reglovh, FUN=function(x) sum(!is.na(x)) ))
	abyti_reglovl<- as.matrix(aggregate(support_12 ~ treatment, data=d2_reglovl, FUN=function(x) sum(!is.na(x)) ))

		a1_reglovh<-abyti_reglovh[1,2]
		a2_reglovh<-abyti_reglovh[2,2]
		a3_reglovh<-abyti_reglovh[3,2]
		a4_reglovh<-abyti_reglovh[4,2]
		a5_reglovh<-abyti_reglovh[5,2]

		a1_reglovl<-abyti_reglovl[1,2]
		a2_reglovl<-abyti_reglovl[2,2]
		a3_reglovl<-abyti_reglovl[3,2]
		a4_reglovl<-abyti_reglovl[4,2]
		a5_reglovl<-abyti_reglovl[5,2]
	
	# Number of draws from the Beta posterior
		n<-5000
	# Draws from the posterior for each treatment
		pi_1_reglovh <- as.matrix(rbeta(n, a1_reglovh+0.5, n1_reglovh-a1_reglovh+0.5))
		pi_2_reglovh <- as.matrix(rbeta(n, a2_reglovh+0.5, n2_reglovh-a2_reglovh+0.5))
		pi_3_reglovh <- as.matrix(rbeta(n, a3_reglovh+0.5, n3_reglovh-a3_reglovh+0.5))
		pi_4_reglovh <- as.matrix(rbeta(n, a4_reglovh+0.5, n4_reglovh-a4_reglovh+0.5))
		pi_5_reglovh <- as.matrix(rbeta(n, a5_reglovh+0.5, n5_reglovh-a5_reglovh+0.5))

		pi_1_reglovl <- as.matrix(rbeta(n, a1_reglovl+0.5, n1_reglovl-a1_reglovl+0.5))
		pi_2_reglovl <- as.matrix(rbeta(n, a2_reglovl+0.5, n2_reglovl-a2_reglovl+0.5))
		pi_3_reglovl <- as.matrix(rbeta(n, a3_reglovl+0.5, n3_reglovl-a3_reglovl+0.5))
		pi_4_reglovl <- as.matrix(rbeta(n, a4_reglovl+0.5, n4_reglovl-a4_reglovl+0.5))
		pi_5_reglovl <- as.matrix(rbeta(n, a5_reglovl+0.5, n5_reglovl-a5_reglovl+0.5))

	# Treatment indicators to make the data for the plot cleaner
		t1 <- as.matrix(rep(1,n, nrow=n, ncol=1))
		t2 <- as.matrix(rep(2,n, nrow=n, ncol=1))
		t3 <- as.matrix(rep(3,n, nrow=n, ncol=1))
		t4 <- as.matrix(rep(4,n, nrow=n, ncol=1))
		t5 <- as.matrix(rep(5,n, nrow=n, ncol=1))

	tis <- as.matrix(rbind(t1,t2,t3,t4,t5))

	pis_reglovh <- as.matrix(rbind(pi_1_reglovh,pi_2_reglovh,pi_3_reglovh,pi_4_reglovh,pi_5_reglovh))
	pis_reglovl <- as.matrix(rbind(pi_1_reglovl,pi_2_reglovl,pi_3_reglovl,pi_4_reglovl,pi_5_reglovl))

	betas_reglovh <- as.data.frame(cbind(tis,pis_reglovh))
	betas_reglovl <- as.data.frame(cbind(tis,pis_reglovl))

	betas_reglovh$treatmentname = factor(betas_reglovh$V1, labels = c("Direct","Interdep.","Network","Null","Plac."))
	betas_reglovl$treatmentname = factor(betas_reglovl$V1, labels = c("Direct","Interdep.","Network","Null","Plac."))

	###
	# Figure 2 (1 of 6 -- top right)
	###
	#pdf("PercSupport_12_Mod_RegLove_BetaEsts.pdf")
		# Approval ratings by treatment group, reglovh
	lineplot.CI(x.factor = betas_reglovh$treatmentname, response = betas_reglovh$V2, type="p",ylab="",xlab="",axes=F, font=2,
		data = betas_reglovh, #main = "Perceive `Too Little' Regulation", #xlab = "Treatment Group", ylab = "% Supporting Regulation (Medium)",
		ylim=c(0,1),
		ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)) )
	#mtext("Percent Supporting More Financial Regulation\n ``too little current regulation''", side=3,cex=1.15, font=2, line=1)
	mtext("``Too little'' current regulation", side=3,cex=2, font=2, line=1)
	mtext("(n=427)", side=3,cex=1.25, font=2, line=0)
		axis(2, las=2)
		mtext("Treatment Group", side=1,cex=1.15, font=1, line=2.5)
	#dev.off()
	###
	# Figure 2 (2 of 6 -- top left)
	###
	#pdf("PercSupport_12_Mod_RegHate_BetaEsts.pdf")
		# Approval ratings by treatment group, reglovh
	lineplot.CI(x.factor = betas_reglovl$treatmentname, response = betas_reglovl$V2, type="p",ylab="",xlab="",axes=F, font=2,
		data = betas_reglovl, #main = "Perceive `Too Much' Regulation", #xlab = "Treatment Group", ylab = "% Supporting Regulation (Medium)",
		ylim=c(0,1),
		ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)) )
		#mtext("Percent Supporting More Financial Regulation\n ``too much current regulation''", side=3,cex=1.15, font=2, line=1)
	mtext("``Too much'' current regulation", side=3,cex=2, font=2, line=1)
	mtext("(n=235)", side=3,cex=1.25, font=2, line=0)
		axis(2, las=2)
		mtext("Treatment Group", side=1,cex=1.15, font=1, line=2.5)
	#dev.off()

	###
	# Fig 2, graphs 3 and 4 (of 6) - PARTISANSHIP
	###
		d_partyDem <- d[ which(d$party_1 == "Democrat"), ]	#Democrat
		d_partyRep <- d[ which(d$party_1 == "Republican"), ]	#Rep
	d3$party_clear <- rep(NA,nrow(d3))
		unique(d3$party_1)
	d3$party_clear[d3$party_1=="Republican"] <- 1
	d3$party_clear[d3$party_1=="Democrat"] <- 3
	###
	# Support 12
	# Treatment pictures for Rep / Dem
	###
	# Summarizes the number of respondents that received each treatment, by h/l
	nbyti_partyRep <- as.matrix(aggregate(support_12 ~ treatment , data=d_partyRep, FUN=function(x) sum(!is.na(x)) ))
	nbyti_partyDem <- as.matrix(aggregate(support_12 ~ treatment , data=d_partyDem, FUN=function(x) sum(!is.na(x)) ))
		n1_partyRep<-nbyti_partyRep[1,2]
		n2_partyRep<-nbyti_partyRep[2,2]
		n3_partyRep<-nbyti_partyRep[3,2]
		n4_partyRep<-nbyti_partyRep[4,2]
		n5_partyRep<-nbyti_partyRep[5,2]

		n1_partyDem<-nbyti_partyDem[1,2]
		n2_partyDem<-nbyti_partyDem[2,2]
		n3_partyDem<-nbyti_partyDem[3,2]
		n4_partyDem<-nbyti_partyDem[4,2]
		n5_partyDem<-nbyti_partyDem[5,2]

	# Repeats that step, for the number of respondents by treatment who approved
	d2_partyRep <- d_partyRep[ which(d_partyRep$support_12 =="1"), ]
	d2_partyDem <- d_partyDem[ which(d_partyDem$support_12 =="1"), ]

	abyti_partyRep<- as.matrix(aggregate(support_12 ~ treatment, data=d2_partyRep, FUN=function(x) sum(!is.na(x)) ))
	abyti_partyDem<- as.matrix(aggregate(support_12 ~ treatment, data=d2_partyDem, FUN=function(x) sum(!is.na(x)) ))

		a1_partyDem<-abyti_partyDem[1,2]
		a2_partyDem<-abyti_partyDem[2,2]
		a3_partyDem<-abyti_partyDem[3,2]
		a4_partyDem<-abyti_partyDem[4,2]
		a5_partyDem<-abyti_partyDem[5,2]

		a1_partyRep<-abyti_partyRep[1,2]
		a2_partyRep<-abyti_partyRep[2,2]
		a3_partyRep<-abyti_partyRep[3,2]
		a4_partyRep<-abyti_partyRep[4,2]
		a5_partyRep<-abyti_partyRep[5,2]

	# Number of draws from the Beta posterior
		n<-5000

	# Draws from the posterior for each treatment
		pi_1_partyRep <- as.matrix(rbeta(n, a1_partyRep+0.5, n1_partyRep-a1_partyRep+0.5))
		pi_2_partyRep <- as.matrix(rbeta(n, a2_partyRep+0.5, n2_partyRep-a2_partyRep+0.5))
		pi_3_partyRep <- as.matrix(rbeta(n, a3_partyRep+0.5, n3_partyRep-a3_partyRep+0.5))
		pi_4_partyRep <- as.matrix(rbeta(n, a4_partyRep+0.5, n4_partyRep-a4_partyRep+0.5))
		pi_5_partyRep <- as.matrix(rbeta(n, a5_partyRep+0.5, n5_partyRep-a5_partyRep+0.5))

		pi_1_partyDem <- as.matrix(rbeta(n, a1_partyDem+0.5, n1_partyDem-a1_partyDem+0.5))
		pi_2_partyDem <- as.matrix(rbeta(n, a2_partyDem+0.5, n2_partyDem-a2_partyDem+0.5))
		pi_3_partyDem <- as.matrix(rbeta(n, a3_partyDem+0.5, n3_partyDem-a3_partyDem+0.5))
		pi_4_partyDem <- as.matrix(rbeta(n, a4_partyDem+0.5, n4_partyDem-a4_partyDem+0.5))
		pi_5_partyDem <- as.matrix(rbeta(n, a5_partyDem+0.5, n5_partyDem-a5_partyDem+0.5))

	# Treatment indicators to make the data for the plot cleaner
		t1 <- as.matrix(rep(1,n, nrow=n, ncol=1))
		t2 <- as.matrix(rep(2,n, nrow=n, ncol=1))
		t3 <- as.matrix(rep(3,n, nrow=n, ncol=1))
		t4 <- as.matrix(rep(4,n, nrow=n, ncol=1))
		t5 <- as.matrix(rep(5,n, nrow=n, ncol=1))

	tis <- as.matrix(rbind(t1,t2,t3,t4,t5))

	pis_partyRep <- as.matrix(rbind(pi_1_partyRep,pi_2_partyRep,pi_3_partyRep,pi_4_partyRep,pi_5_partyRep))
	pis_partyDem <- as.matrix(rbind(pi_1_partyDem,pi_2_partyDem,pi_3_partyDem,pi_4_partyDem,pi_5_partyDem))

	betas_partyRep <- as.data.frame(cbind(tis,pis_partyRep))
	betas_partyDem <- as.data.frame(cbind(tis,pis_partyDem))

	betas_partyRep$treatmentname = factor(betas_partyRep$V1, labels = c("Direct","Interdep.","Network","Null","Plac."))
	betas_partyDem$treatmentname = factor(betas_partyDem$V1, labels = c("Direct","Interdep.","Network","Null","Plac."))

	###
	# Figure 2 (3 of 6 -- middle left)
	###
	#pdf("PercSupport_12_Mod_PartyRep_BetaEsts.pdf")
		# Approval ratings by treatment group
	lineplot.CI(x.factor = betas_partyRep$treatmentname, response = betas_partyRep$V2, type="p",ylab="",xlab="",axes=F, font=2,
		data = betas_partyRep, 
		ylim=c(0,1),
		ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)) )
	#mtext("Percent Supporting More Financial Regulation\n Republican Partisan Identity", side=3,cex=1.15, font=2, line=1)
	mtext("``Republican''", side=3,cex=2, font=2, line=1)
	mtext("(n=180)", side=3,cex=1.25, font=2, line=0)
		axis(2, las=2)
		mtext("Treatment Group", side=1,cex=1.15, font=1, line=2.5)
	#dev.off()
	###
	# Figure 2 (4 of 6 -- middle right)
	###
	#pdf("PercSupport_12_Mod_PartyDem_BetaEsts.pdf")
		# Approval ratings by treatment group
	lineplot.CI(x.factor = betas_partyDem$treatmentname, response = betas_partyDem$V2, type="p",ylab="",xlab="",axes=F, font=2,
		data = betas_partyDem, 
		ylim=c(0,1),
		ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)) )
	#mtext("Percent Supporting More Financial Regulation\n Democratic Partisan Identity", side=3,cex=1.15, font=2, line=1)
	mtext("``Democrat''", side=3,cex=2, font=2, line=1)
	mtext("(n=465)", side=3,cex=1.25, font=2, line=0)
		axis(2, las=2)
		mtext("Treatment Group", side=1,cex=1.15, font=1, line=2.5)
	#dev.off()

	###
	# Fig 2, graphs 5 and 6 (of 6) - IDEOLOGY
	###
		d_ideolLib <- d[ which(d$libcon == "Extremely liberal"|
						d$libcon == "Liberal" |
						d$libcon == "Slightly liberal"
					), ]	#Liberal
		d_ideolCons<- d[ which(d$libcon == "Extremely conservative"|
						d$libcon == "Conservative" |
						d$libcon == "Slightly conservative"
					), ]	#Conservative
	# Summarizes the number of respondents that received each treatment, by h/l
	nbyti_ideolCons <- as.matrix(aggregate(support_12 ~ treatment , data=d_ideolCons, FUN=function(x) sum(!is.na(x)) ))
	nbyti_ideolLib <- as.matrix(aggregate(support_12 ~ treatment , data=d_ideolLib, FUN=function(x) sum(!is.na(x)) ))

		n1_ideolCons<-nbyti_ideolCons[1,2]
		n2_ideolCons<-nbyti_ideolCons[2,2]
		n3_ideolCons<-nbyti_ideolCons[3,2]
		n4_ideolCons<-nbyti_ideolCons[4,2]
		n5_ideolCons<-nbyti_ideolCons[5,2]

		n1_ideolLib<-nbyti_ideolLib[1,2]
		n2_ideolLib<-nbyti_ideolLib[2,2]
		n3_ideolLib<-nbyti_ideolLib[3,2]
		n4_ideolLib<-nbyti_ideolLib[4,2]
		n5_ideolLib<-nbyti_ideolLib[5,2]

	# Repeats that step, for the number of respondents by treatment who approved
	d2_ideolCons <- d_ideolCons[ which(d_ideolCons$support_12 =="1"), ]
	d2_ideolLib <- d_ideolLib[ which(d_ideolLib$support_12 =="1"), ]

	abyti_ideolCons<- as.matrix(aggregate(support_12 ~ treatment, data=d2_ideolCons, FUN=function(x) sum(!is.na(x)) ))
	abyti_ideolLib<- as.matrix(aggregate(support_12 ~ treatment, data=d2_ideolLib, FUN=function(x) sum(!is.na(x)) ))

		a1_ideolLib<-abyti_ideolLib[1,2]
		a2_ideolLib<-abyti_ideolLib[2,2]
		a3_ideolLib<-abyti_ideolLib[3,2]
		a4_ideolLib<-abyti_ideolLib[4,2]
		a5_ideolLib<-abyti_ideolLib[5,2]

		a1_ideolCons<-abyti_ideolCons[1,2]
		a2_ideolCons<-abyti_ideolCons[2,2]
		a3_ideolCons<-abyti_ideolCons[3,2]
		a4_ideolCons<-abyti_ideolCons[4,2]
		a5_ideolCons<-abyti_ideolCons[5,2]
	
	# Number of draws from the Beta posterior
		n<-5000

	# Draws from the posterior for each treatment
		pi_1_ideolCons <- as.matrix(rbeta(n, a1_ideolCons+0.5, n1_ideolCons-a1_ideolCons+0.5))
		pi_2_ideolCons <- as.matrix(rbeta(n, a2_ideolCons+0.5, n2_ideolCons-a2_ideolCons+0.5))
		pi_3_ideolCons <- as.matrix(rbeta(n, a3_ideolCons+0.5, n3_ideolCons-a3_ideolCons+0.5))
		pi_4_ideolCons <- as.matrix(rbeta(n, a4_ideolCons+0.5, n4_ideolCons-a4_ideolCons+0.5))
		pi_5_ideolCons <- as.matrix(rbeta(n, a5_ideolCons+0.5, n5_ideolCons-a5_ideolCons+0.5))

		pi_1_ideolLib <- as.matrix(rbeta(n, a1_ideolLib+0.5, n1_ideolLib-a1_ideolLib+0.5))
		pi_2_ideolLib <- as.matrix(rbeta(n, a2_ideolLib+0.5, n2_ideolLib-a2_ideolLib+0.5))
		pi_3_ideolLib <- as.matrix(rbeta(n, a3_ideolLib+0.5, n3_ideolLib-a3_ideolLib+0.5))
		pi_4_ideolLib <- as.matrix(rbeta(n, a4_ideolLib+0.5, n4_ideolLib-a4_ideolLib+0.5))
		pi_5_ideolLib <- as.matrix(rbeta(n, a5_ideolLib+0.5, n5_ideolLib-a5_ideolLib+0.5))

	# Treatment indicators to make the data for the plot cleaner
		t1 <- as.matrix(rep(1,n, nrow=n, ncol=1))
		t2 <- as.matrix(rep(2,n, nrow=n, ncol=1))
		t3 <- as.matrix(rep(3,n, nrow=n, ncol=1))
		t4 <- as.matrix(rep(4,n, nrow=n, ncol=1))
		t5 <- as.matrix(rep(5,n, nrow=n, ncol=1))

	tis <- as.matrix(rbind(t1,t2,t3,t4,t5))

	pis_ideolCons <- as.matrix(rbind(pi_1_ideolCons,pi_2_ideolCons,pi_3_ideolCons,pi_4_ideolCons,pi_5_ideolCons))
	pis_ideolLib <- as.matrix(rbind(pi_1_ideolLib,pi_2_ideolLib,pi_3_ideolLib,pi_4_ideolLib,pi_5_ideolLib))

	betas_ideolCons <- as.data.frame(cbind(tis,pis_ideolCons))
	betas_ideolLib <- as.data.frame(cbind(tis,pis_ideolLib))

	betas_ideolCons$treatmentname = factor(betas_ideolCons$V1, labels = c("Direct","Interdep.","Network","Null","Plac."))
	betas_ideolLib$treatmentname = factor(betas_ideolLib$V1, labels = c("Direct","Interdep.","Network","Null","Plac."))
	
	###
	# Figure 2 (5 of 6 -- bottom right)
	###
	#pdf("PercSupport_12_Mod_ideolLib_BetaEsts.pdf")
		# Approval ratings by treatment group
	lineplot.CI(x.factor = betas_ideolLib$treatmentname, response = betas_ideolLib$V2, type="p",ylab="",xlab="",axes=F, font=2,
		data = betas_ideolLib, 
		ylim=c(0,1),
		ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)) )
	#mtext("Percent Supporting More Financial Regulation\n Liberal Ideology", side=3,cex=1.15, font=2, line=1)
	mtext("``Liberal''", side=3,cex=2, font=2, line=1)
	mtext("(n=580)", side=3,cex=1.25, font=2, line=0)
		axis(2, las=2)
		mtext("Treatment Group", side=1,cex=1.15, font=1, line=2.5)
	#dev.off()

	###
	# Figure 2 (6 of 6 -- bottom left)
	###
	#pdf("PercSupport_12_Mod_ideolCons_BetaEsts.pdf")
		# Approval ratings by treatment group
	lineplot.CI(x.factor = betas_ideolCons$treatmentname, response = betas_ideolCons$V2, type="p",ylab="",xlab="",axes=F, font=2,
		data = betas_ideolCons, 
		ylim=c(0,1),
		ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)) )
	#mtext("Percent Supporting More Financial Regulation\n Conservative Ideology", side=3,cex=1.15, font=2, line=1)
	mtext("``Conservative''", side=3,cex=2, font=2, line=1)
	mtext("(n=249)", side=3,cex=1.25, font=2, line=0)
		axis(2, las=2)
		mtext("Treatment Group", side=1,cex=1.15, font=1, line=2.5)
	#dev.off()





###
# Table 2
###
	#######
	#LOGIT analysis -- Table 2
	#######
	d <- read.csv("ChaudoinWilf_FPA_data.csv")
	nrow(d3<-subset(d, v10==1 & consent_tx=="Yes" & is.na(support_12)==F & minutes>=8 & minutes<=32)) #n = 1049
	data <- d3
	n <- length(data$support_12)
	y_12 <- data$support_12
	y_123 <- data$support_123
	y_1234 <- data$support_1234

	treat_oep <- data$treat_oep 
	treat_nia <- data$treat_nia 
	treat_nw <- data$treat_nw 
	treat_plac <- data$treat_plac 

	#ADD OTHER VARIABLES. . . 
		data$education_7ptscale <-rep(NA,nrow(data))
		data$education_7ptscale[data$education=="9th - 12th grade (no diploma or GED)"] <-1
			sum(data$education_7ptscale==1,na.rm=T)		#3
		data$education_7ptscale[data$education=="High School graduate (diploma or GED)"] <-2
			sum(data$education_7ptscale==2,na.rm=T)		#101
		data$education_7ptscale[data$education=="Some college, no degree"] <-3
			sum(data$education_7ptscale==3,na.rm=T)		#271
		data$education_7ptscale[data$education=="Associates degree"] <-4
			sum(data$education_7ptscale==4,na.rm=T)		#108
		data$education_7ptscale[data$education=="Bachelors degree"] <-5
			sum(data$education_7ptscale==5,na.rm=T)		#359
		data$education_7ptscale[data$education=="Masters degree"] <-6
			sum(data$education_7ptscale==6,na.rm=T)		# 82
		data$education_7ptscale[data$education=="Professional or Doctorate degree"] <-7
			sum(data$education_7ptscale==7,na.rm=T)		# 22
	edu7pt <- data$education_7ptscale
	edu <- data$educ_ba	#BINARY education

	white <- data$white; unique(data$white)
	female <- data$female_ind; unique(data$female_ind)
	employed <- data$employed; unique(data$employed)
	incomegt50k <- data$income_gt50k; unique(data$income_gt50k)
		income_10pt <- data$income_num; unique(data$income_num)
	pk <- data$pk_sum; unique(data$pk_sum)

	dem <- data$democrat; unique(data$democrat)
	ideology <- data$democrat; unique(data$democrat)
		data$ideology_7ptscale <-rep(NA,nrow(data))
		data$ideology_7ptscale[data$libcons=="Extremely conservative"] <-1
			sum(data$ideology_7ptscale==1,na.rm=T)		#30
		data$ideology_7ptscale[data$libcons=="Conservative"] <-2
			sum(data$ideology_7ptscale==2,na.rm=T)		#105
		data$ideology_7ptscale[data$libcons=="Slightly conservative"] <-3
			sum(data$ideology_7ptscale==3,na.rm=T)		#85
		data$ideology_7ptscale[data$libcons=="Moderate, middle of the road"] <-4
			sum(data$ideology_7ptscale==4,na.rm=T)		#196
		data$ideology_7ptscale[data$libcons=="Slightly liberal"] <-5
			sum(data$ideology_7ptscale==5,na.rm=T)		#147
		data$ideology_7ptscale[data$libcons=="Liberal"] <-6
			sum(data$ideology_7ptscale==6,na.rm=T)		#267
		data$ideology_7ptscale[data$libcons=="Extremely liberal"] <-7
			sum(data$ideology_7ptscale==7,na.rm=T)		#115
	ideology7pt <- data$ideology_7ptscale
		data$reglove3pt <- rep(NA, nrow(data))
		data$reglove3pt[data$reglove1=="Too much"]  <- 1; sum(data$reglove3pt==1,na.rm=T)
		data$reglove3pt[data$reglove1=="About right"]  <- 2; sum(data$reglove3pt==2,na.rm=T)
		data$reglove3pt[data$reglove1=="Too little"]  <- 3; sum(data$reglove3pt==3,na.rm=T)
	reglove3pt <- data$reglove3pt

	fr3pt <- data$frsum_red
	nonfr_ind <- data$frsum_red
		data$fr_123_ind <- rep(NA, nrow(data))
		data$fr_123_ind[data$frsum_red>0]  <- 1; sum(data$fr_123_ind==1,na.rm=T)
		data$fr_123_ind[data$frsum_red==0]  <- 0; sum(data$fr_123_ind==0,na.rm=T)
	fr_123_ind <- data$fr_123_ind
		data$fr_0_ind <- rep(NA, nrow(data))
		data$fr_0_ind[data$frsum_red>0]  <- 0; sum(data$fr_0_ind==0,na.rm=T)
		data$fr_0_ind[data$frsum_red==0]  <- 1; sum(data$fr_0_ind==1,na.rm=T)
	fr_0_ind <- data$fr_0_ind
		data$fr_23_ind <- rep(NA, nrow(data))
		data$fr_23_ind[data$frsum_red>1]  <- 1; sum(data$fr_23_ind==1,na.rm=T)
		data$fr_23_ind[data$frsum_red<2]  <- 0; sum(data$fr_23_ind==0,na.rm=T)
	fr_23_ind <- data$fr_23_ind
	ethno_val <- data$ethno_diff
		(q2val<-as.numeric(summary(data$ethno_diff)[3]))
			summary(data$ethno_diff)	#min-max = -6.25 to 3.361; mean=-1.187, 1q/2q/3q = -1.583, -1.22, -1.00
		data$ansamp_ethnoh_ind <- ifelse(is.na(data$ethno_diff)==F,ifelse(data$ethno_diff>q2val,1,0),NA)
	ethno_high_ind <- data$ansamp_ethnoh_ind

	## LOGITS 
	#Table 2, Model (1)
	logit1 <-logit <- glm(y_12 ~ treat_oep + treat_nia + treat_nw + treat_plac,
			data = data, family = binomial)
		summary(logit)
	#Table 2, Model (2)
	logit2 <-logit <- glm(y_12 ~ treat_oep + treat_nia + treat_nw + treat_plac
			+ edu7pt + white + female + employed + incomegt50k + pk, 
			data = data, family = binomial)
		summary(logit)
		(nrow(sub <- subset(data, is.na(y_12)==F 
						& is.na(treat_oep)==F
						& is.na(treat_nia)==F
						& is.na(treat_nw)==F
						& is.na(treat_plac)==F
						& is.na(edu7pt)==F
						& is.na(white)==F
						& is.na(female)==F
						& is.na(employed)==F
						& is.na(income_gt50k)==F
						& is.na(pk_sum)==F))) #n=1045
	#Table 2, Model (3)
	logit3 <-logit <- glm(y_12 ~ treat_oep + treat_nia + treat_nw + treat_plac
			+ edu7pt + white + female + employed + incomegt50k + pk 
			+ reglove3pt, 
			data = data, family = binomial)
		summary(logit)
		(nrow(sub <- subset(data, is.na(y_12)==F 
						& is.na(treat_oep)==F
						& is.na(treat_nia)==F
						& is.na(treat_nw)==F
						& is.na(treat_plac)==F
						& is.na(edu7pt)==F
						& is.na(white)==F
						& is.na(female)==F
						& is.na(employed)==F
						& is.na(incomegt50k)==F
						& is.na(pk)==F
						& is.na(reglove3pt)==F)))	#n=1044
	#Table 2, Model (4)
	logit4 <-logit <- glm(y_12 ~ treat_oep + treat_nia + treat_nw + treat_plac
			+ edu7pt + white + female + employed + incomegt50k + pk + democrat, 
			data = data, family = binomial)
		summary(logit)

		(nrow(sub <- subset(data, is.na(y_12)==F 
						& is.na(treat_oep)==F
						& is.na(treat_nia)==F
						& is.na(treat_nw)==F
						& is.na(treat_plac)==F
						& is.na(edu7pt)==F
						& is.na(white)==F
						& is.na(female)==F
						& is.na(employed)==F
						& is.na(incomegt50k)==F
						& is.na(pk)==F
						& is.na(democrat)==F)))	#n=1044
						#& is.na(reglove3pt)==F)))	#n=1044
	#Table 2, Model (5)
	logit5 <-logit <- glm(y_12 ~ treat_oep + treat_nia + treat_nw + treat_plac
			+ edu7pt + white + female + employed + incomegt50k + pk 
			+ ideology7pt, 
			data = data, family = binomial)
		summary(logit)
	(nrow(sub <- subset(data, is.na(y_12)==F 
						& is.na(treat_oep)==F
						& is.na(treat_nia)==F
						& is.na(treat_nw)==F
						& is.na(treat_plac)==F
						& is.na(edu7pt)==F
						& is.na(white)==F
						& is.na(female)==F
						& is.na(employed)==F
						& is.na(incomegt50k)==F
						& is.na(pk)==F
						& is.na(ideology_7ptscale)==F)))	#n=1044
						#& is.na(democrat)==F)))	#n=1044
						#& is.na(reglove3pt)==F)))	#n=1044
	#Table 2, Model (6)
	logit6 <-logit <- glm(y_12 ~ treat_oep + treat_nia + treat_nw + treat_plac
			+ edu7pt + white + female + employed + incomegt50k + pk + ideology7pt
			+ fr3pt + ethno_val, 
			#+ fr3pt + ethno_high_ind, 
			data = data, family = binomial)
		summary(logit)
	(nrow(sub <- subset(data, is.na(y_12)==F 
						& is.na(treat_oep)==F
						& is.na(treat_nia)==F
						& is.na(treat_nw)==F
						& is.na(treat_plac)==F
						& is.na(edu7pt)==F
						& is.na(white)==F
						& is.na(female)==F
						& is.na(employed)==F
						& is.na(incomegt50k)==F
						& is.na(pk)==F
						& is.na(ideology_7ptscale)==F
						& is.na(frsum_red)==F
						& is.na(ethno_val)==F)))	# 900

	install.packages("stargazer")
	library(stargazer)
	stargazer(logit1, logit2, logit3, logit4, logit5, logit6, style="qje",
			align=T, ci=F,  	dep.var.labels.include = FALSE,
			no.space=TRUE)


############################
# Table B.1
############################
	#########################
	# Moderator Descriptive Statistics
	#########################

	#1049 respondents in the sample
	d <- read.csv("ChaudoinWilf_FPA_data.csv")
	nrow(d3<-subset(d, v10==1 & consent_tx=="Yes" & is.na(support_12)==F & minutes>=8 & minutes<=32)) #n = 1049


	############################
	#Folk realism values
	############################
	#Pooled
	(frval <- sum(is.na(d3$frsum_red)==F))	#1040 respondents have FR moderator values
		#34.9% never chose the realist option
			sum(d3$frsum_red==0,na.rm=T)/frval
		#65.1% chose the realist option 1+ time
			sum(d3$frsum_red==1 | d3$frsum_red==2 | d3$frsum_red==3,na.rm=T)/frval
	#Direct
		(denom <- nrow(subset(d3, is.na(frsum_red)==F & treat_oep==1)))
		#38.3% oep-treated never chose the realist option
			nrow(subset(d3,frsum_red==0 & treat_oep==1,na.rm=T))/denom
		#61.7% oep-treated never chose the realist option
			nrow(subset(d3,(frsum_red==1|frsum_red==2|frsum_red==3) & treat_oep==1,na.rm=T))/denom
	#Interdependent
		(denom <- nrow(subset(d3, is.na(frsum_red)==F & treat_nia==1)))
		#32.5% interdep.-treated never chose the realist option
			nrow(subset(d3,frsum_red==0 & treat_nia==1,na.rm=T))/denom
		#67.5% interdep.-treated never chose the realist option
			nrow(subset(d3,(frsum_red==1|frsum_red==2|frsum_red==3) & treat_nia==1,na.rm=T))/denom
	#Network
		(denom <- nrow(subset(d3, is.na(frsum_red)==F & treat_nw==1)))
		#32.2% network-treated never chose the realist option
			nrow(subset(d3,frsum_red==0 & treat_nw==1,na.rm=T))/denom
		#67.8% network-treated never chose the realist option
			nrow(subset(d3,(frsum_red==1|frsum_red==2|frsum_red==3) & treat_nw==1,na.rm=T))/denom
	#Null
		(denom <- nrow(subset(d3, is.na(frsum_red)==F & treat_null==1)))
		#32.2% null arg recipients never chose the realist option
			nrow(subset(d3,frsum_red==0 & treat_null==1,na.rm=T))/denom
		#67.8% null arg recipients never chose the realist option
			nrow(subset(d3,(frsum_red==1|frsum_red==2|frsum_red==3) & treat_null==1,na.rm=T))/denom
	#Placebo
		(denom <- nrow(subset(d3, is.na(frsum_red)==F & treat_plac==1)))
		#39.3% placebo recipients never chose the realist option
			nrow(subset(d3,frsum_red==0 & treat_plac==1,na.rm=T))/denom
		#60.7% placebo recipients never chose the realist option
			nrow(subset(d3,(frsum_red==1|frsum_red==2|frsum_red==3) & treat_plac==1,na.rm=T))/denom

	############################
	#Ethnocentrism
	############################

	mv <- medianVal <- summary(d3$ethno_diff)[3]
	#Pooled
		(denom <- nrow(subset(d3, is.na(ethno_diff)==F)))
		#49.1% at and above median
			(e11<-sum(d3$ethno_diff>=mv,na.rm=T)/denom)
		#50.9% below median
			(e12<-sum(d3$ethno_diff<mv,na.rm=T)/denom)
	#Direct-treatment
		(denom <- nrow(subset(d3, is.na(ethno_diff)==F & treat_oep==1)))
		#52.4% at and above median
			(e21<-sum(d3$ethno_diff>=mv & d3$treat_oep==1,na.rm=T)/denom)
		#47.6% below median
			(e22<-sum(d3$ethno_diff<mv& d3$treat_oep==1,na.rm=T)/denom)
	#Interdep. Treatment
		(denom <- nrow(subset(d3, is.na(ethno_diff)==F & treat_nia==1)))
		#51.4% at and above median
			(e31<-sum(d3$ethno_diff>=mv & d3$treat_nia==1,na.rm=T)/denom)
		#48.6% below median
			(e32<-sum(d3$ethno_diff<mv& d3$treat_nia==1,na.rm=T)/denom)
	#Network Treatment
		(denom <- nrow(subset(d3, is.na(ethno_diff)==F & treat_nw==1)))
		#48.4% at and above median
			(e41<-sum(d3$ethno_diff>=mv & d3$treat_nw==1,na.rm=T)/denom)
		#51.6% below median
			(e42<-sum(d3$ethno_diff<mv& d3$treat_nw==1,na.rm=T)/denom)
	#Null argument
		(denom <- nrow(subset(d3, is.na(ethno_diff)==F & treat_null==1)))
		#46.6% at and above median
			(e51<-sum(d3$ethno_diff>=mv & d3$treat_null==1,na.rm=T)/denom)
		#53.4% below median
			(e52<-sum(d3$ethno_diff<mv& d3$treat_null==1,na.rm=T)/denom)
	#Placebo argument
		(denom <- nrow(subset(d3, is.na(ethno_diff)==F & treat_plac==1)))
		#46.5% at and above median
			(e61<-sum(d3$ethno_diff>=mv & d3$treat_plac==1,na.rm=T)/denom)
		#53.5% below median
			(e62<-sum(d3$ethno_diff<mv& d3$treat_plac==1,na.rm=T)/denom)

	round(c(e11,e21,e31,e41,e51,e61)*100,1)#Ethnocentrist percent per treatment
	round(c(e12,e22,e32,e42,e52,e62)*100,1)#Non-Ethnocentrist percent per treatment

############################
# Table B.2
############################
	d <- read.csv("ChaudoinWilf_FPA_data.csv")
	nrow(d3<-subset(d, v10==1 & consent_tx=="Yes" & is.na(support_12)==F & minutes>=8 & minutes<=32)) #n = 1049

	#Pooled
	(AllRespondents <- nrow(subset(d,  treat_oep==1|treat_nia==1|treat_nw==1|treat_null==1|treat_plac==1)))
	(InSampleRespondents <- nrow(subset(d3, v10==1)))
	(Support <- nrow(subset(d3, support_12==1))/InSampleRespondents)
	(SupportModerate <- nrow(subset(d3, support_123==1))/InSampleRespondents)
	(SupportBroad <- nrow(subset(d3,  support_1234==1))/InSampleRespondents)

	#Direct
	(AllRespondents <- nrow(subset(d,  treat_oep==1)))
	(InSampleRespondents <- nrow(subset(d3,  treat_oep==1)))
	(Support <- nrow(subset(d3,  treat_oep==1 & support_12==1))/InSampleRespondents)
	(SupportModerate <- nrow(subset(d3,  treat_oep==1 & support_123==1))/InSampleRespondents)
	(SupportBroad <- nrow(subset(d3,  treat_oep==1 & support_1234==1))/InSampleRespondents)

	#Interdependence
	(AllRespondents <- nrow(subset(d,  treat_nia==1)))
	(InSampleRespondents <- nrow(subset(d3,  treat_nia==1)))
	(Support <- nrow(subset(d3,  treat_nia==1 & support_12==1))/InSampleRespondents)
	(SupportModerate <- nrow(subset(d3,  treat_nia==1 & support_123==1))/InSampleRespondents)
	(SupportBroad <- nrow(subset(d3,  treat_nia==1 & support_1234==1))/InSampleRespondents)

	#Network
	(AllRespondents <- nrow(subset(d,  treat_nw==1)))
	(InSampleRespondents <- nrow(subset(d3,  treat_nw==1)))
	(Support <- nrow(subset(d3,  treat_nw==1 & support_12==1))/InSampleRespondents)
	(SupportModerate <- nrow(subset(d3,  treat_nw==1 & support_123==1))/InSampleRespondents)
	(SupportBroad <- nrow(subset(d3,  treat_nw==1 & support_1234==1))/InSampleRespondents)

	#Null
	(AllRespondents <- nrow(subset(d,  treat_null==1)))
	(InSampleRespondents <- nrow(subset(d3,  treat_null==1)))
	(Support <- nrow(subset(d3,  treat_null==1 & support_12==1))/InSampleRespondents)
	(SupportModerate <- nrow(subset(d3,  treat_null==1 & support_123==1))/InSampleRespondents)
	(SupportBroad <- nrow(subset(d3,  treat_null==1 & support_1234==1))/InSampleRespondents)

	#Placebo
	(AllRespondents <- nrow(subset(d,  treat_plac==1)))
	(InSampleRespondents <- nrow(subset(d3,  treat_plac==1)))
	(Support <- nrow(subset(d3,  treat_plac==1 & support_12==1))/InSampleRespondents)
	(SupportModerate <- nrow(subset(d3,  treat_plac==1 & support_123==1))/InSampleRespondents)
	(SupportBroad <- nrow(subset(d3,  treat_plac==1 & support_1234==1))/InSampleRespondents)


####
#Table B.3 - Correlation Matrix for Logits
####
	d <- read.csv("ChaudoinWilf_FPA_data.csv")
	nrow(d3<-subset(d, v10==1 & consent_tx=="Yes" & is.na(support_12)==F & minutes>=8 & minutes<=32)) #n = 1049
	data <- d3
	n <- length(data$support_12)
	y_12 <- data$support_12
	y_123 <- data$support_123
	y_1234 <- data$support_1234

	treat_oep <- data$treat_oep 
	treat_nia <- data$treat_nia 
	treat_nw <- data$treat_nw 
	treat_plac <- data$treat_plac 

	#ADD OTHER VARIABLES. . . 
		data$education_7ptscale <-rep(NA,nrow(data))
		data$education_7ptscale[data$education=="9th - 12th grade (no diploma or GED)"] <-1
			sum(data$education_7ptscale==1,na.rm=T)		#3
		data$education_7ptscale[data$education=="High School graduate (diploma or GED)"] <-2
			sum(data$education_7ptscale==2,na.rm=T)		#101
		data$education_7ptscale[data$education=="Some college, no degree"] <-3
			sum(data$education_7ptscale==3,na.rm=T)		#271
		data$education_7ptscale[data$education=="Associates degree"] <-4
			sum(data$education_7ptscale==4,na.rm=T)		#108
		data$education_7ptscale[data$education=="Bachelors degree"] <-5
			sum(data$education_7ptscale==5,na.rm=T)		#359
		data$education_7ptscale[data$education=="Masters degree"] <-6
			sum(data$education_7ptscale==6,na.rm=T)		# 82
		data$education_7ptscale[data$education=="Professional or Doctorate degree"] <-7
			sum(data$education_7ptscale==7,na.rm=T)		# 22
	edu7pt <- data$education_7ptscale
	edu <- data$educ_ba	#BINARY education

	white <- data$white; unique(data$white)
	female <- data$female_ind; unique(data$female_ind)
	employed <- data$employed; unique(data$employed)
	incomegt50k <- data$income_gt50k; unique(data$income_gt50k)
		income_10pt <- data$income_num; unique(data$income_num)
	pk <- data$pk_sum; unique(data$pk_sum)

	dem <- data$democrat; unique(data$democrat)
	ideology <- data$democrat; unique(data$democrat)
		data$ideology_7ptscale <-rep(NA,nrow(data))
		data$ideology_7ptscale[data$libcons=="Extremely conservative"] <-1
			sum(data$ideology_7ptscale==1,na.rm=T)		#30
		data$ideology_7ptscale[data$libcons=="Conservative"] <-2
			sum(data$ideology_7ptscale==2,na.rm=T)		#105
		data$ideology_7ptscale[data$libcons=="Slightly conservative"] <-3
			sum(data$ideology_7ptscale==3,na.rm=T)		#85
		data$ideology_7ptscale[data$libcons=="Moderate, middle of the road"] <-4
			sum(data$ideology_7ptscale==4,na.rm=T)		#196
		data$ideology_7ptscale[data$libcons=="Slightly liberal"] <-5
			sum(data$ideology_7ptscale==5,na.rm=T)		#147
		data$ideology_7ptscale[data$libcons=="Liberal"] <-6
			sum(data$ideology_7ptscale==6,na.rm=T)		#267
		data$ideology_7ptscale[data$libcons=="Extremely liberal"] <-7
			sum(data$ideology_7ptscale==7,na.rm=T)		#115
	ideology7pt <- data$ideology_7ptscale
		data$reglove3pt <- rep(NA, nrow(data))
		data$reglove3pt[data$reglove1=="Too much"]  <- 1; sum(data$reglove3pt==1,na.rm=T)
		data$reglove3pt[data$reglove1=="About right"]  <- 2; sum(data$reglove3pt==2,na.rm=T)
		data$reglove3pt[data$reglove1=="Too little"]  <- 3; sum(data$reglove3pt==3,na.rm=T)
	reglove3pt <- data$reglove3pt

	fr3pt <- data$frsum_red
	nonfr_ind <- data$frsum_red
		data$fr_123_ind <- rep(NA, nrow(data))
		data$fr_123_ind[data$frsum_red>0]  <- 1; sum(data$fr_123_ind==1,na.rm=T)
		data$fr_123_ind[data$frsum_red==0]  <- 0; sum(data$fr_123_ind==0,na.rm=T)
	fr_123_ind <- data$fr_123_ind
		data$fr_0_ind <- rep(NA, nrow(data))
		data$fr_0_ind[data$frsum_red>0]  <- 0; sum(data$fr_0_ind==0,na.rm=T)
		data$fr_0_ind[data$frsum_red==0]  <- 1; sum(data$fr_0_ind==1,na.rm=T)
	fr_0_ind <- data$fr_0_ind
		data$fr_23_ind <- rep(NA, nrow(data))
		data$fr_23_ind[data$frsum_red>1]  <- 1; sum(data$fr_23_ind==1,na.rm=T)
		data$fr_23_ind[data$frsum_red<2]  <- 0; sum(data$fr_23_ind==0,na.rm=T)
	fr_23_ind <- data$fr_23_ind
	ethno_val <- data$ethno_diff
		(q2val<-as.numeric(summary(data$ethno_diff)[3]))
			summary(data$ethno_diff)	#min-max = -6.25 to 3.361; mean=-1.187, 1q/2q/3q = -1.583, -1.22, -1.00
		data$ansamp_ethnoh_ind <- ifelse(is.na(data$ethno_diff)==F,ifelse(data$ethno_diff>q2val,1,0),NA)
	ethno_high_ind <- data$ansamp_ethnoh_ind

	#CORRELATION MATRIX FOR LOGITS
				(nrow(sub <- subset(data, is.na(y_12)==F 
						& is.na(treat_oep)==F
						& is.na(treat_nia)==F
						& is.na(treat_nw)==F
						& is.na(treat_plac)==F
						& is.na(edu7pt)==F
						& is.na(white)==F
						& is.na(female)==F
						& is.na(employed)==F
						& is.na(incomegt50k)==F
						& is.na(pk)==F
						& is.na(ideology_7ptscale)==F
						& is.na(frsum_red)==F
						& is.na(ethno_val)==F	
						& is.na(democrat)==F	
						& is.na(reglove3pt)==F)))	#n=899

	cor(sub[,c("reglove3pt","democrat","ideology_7ptscale","frsum_red","ethno_diff")])


##################
#Table B.4
##################

	d <- read.csv("ChaudoinWilf_FPA_data.csv")
	nrow(d3<-subset(d, v10==1 & consent_tx=="Yes" & is.na(support_12)==F & minutes>=8 & minutes<=32)) #n = 1049

	###############################
	#t-test function code
	t.test2 <- function(m1,m2,s1,s2,n1,n2,m0=0,equal.variance=FALSE)
	{
	    if( equal.variance==FALSE ) 
	    {
        se <- sqrt( (s1^2/n1) + (s2^2/n2) )
        # welch-satterthwaite df
        df <- ( (s1^2/n1 + s2^2/n2)^2 )/( (s1^2/n1)^2/(n1-1) + (s2^2/n2)^2/(n2-1) )
	    } else
    	{
        # pooled standard deviation, scaled by the sample sizes
        se <- sqrt( (1/n1 + 1/n2) * ((n1-1)*s1^2 + (n2-1)*s2^2)/(n1+n2-2) ) 
        df <- n1+n2-2
	    }      
	    t <- (m1-m2-m0)/se 
	    dat <- c(m1-m2, se, t, 2*pt(-abs(t),df))    
	    names(dat) <- c("Difference of means", "Std Error", "t", "p-value")
	    return(dat) 
	}
	###################################

	#DV: Support (Baseline) == 1 or 2
	#Null
	x<-subset(d3, treat_null==1 & is.na(support_12)==F)$support_12
		length(x)	#N
		sum(x==1)/(sum(x==1)+sum(x==0))	#Proportion Support

	#Direct
	y<-subset(d3, treat_oep==1 & is.na(support_12)==F)$support_12
		length(y)	#N
		sum(y==1)/(sum(y==1)+sum(y==0))	#Proportion Support
		t.test2( mean(y), mean(x), sd(y), sd(x), length(y), length(x))	#Diff; SE; T-Stat; p-val

	#Interdep.
	y<-subset(d3, treat_nia==1 & is.na(support_12)==F)$support_12
		length(y)	#N
		sum(y==1)/(sum(y==1)+sum(y==0))	#Proportion Support
		t.test2( mean(y), mean(x), sd(y), sd(x), length(y), length(x))	#Diff; SE; T-Stat; p-val
	#Network
	y<-subset(d3, treat_nw==1 & is.na(support_12)==F)$support_12
		length(y)	#N
		sum(y==1)/(sum(y==1)+sum(y==0))	#Proportion Support
		t.test2( mean(y), mean(x), sd(y), sd(x), length(y), length(x))	#Diff; SE; T-Stat; p-val
	#Placebo
	y<-subset(d3, treat_plac==1 & is.na(support_12)==F)$support_12
		length(y)	#N
		sum(y==1)/(sum(y==1)+sum(y==0))	#Proportion Support
		t.test2( mean(y), mean(x), sd(y), sd(x), length(y), length(x))	#Diff; SE; T-Stat; p-val

	#DV: Support (Moderate Defn) == 1 or 2 or 3
	#Null
	x<-subset(d3, treat_null==1 & is.na(support_123)==F)$support_123
		length(x)	#N
		sum(x==1)/(sum(x==1)+sum(x==0))	#Proportion Support
	#Direct
	y<-subset(d3, treat_oep==1 & is.na(support_123)==F)$support_123
		length(y)	#N
		sum(y==1)/(sum(y==1)+sum(y==0))	#Proportion Support
		t.test2( mean(y), mean(x), sd(y), sd(x), length(y), length(x))	#Diff; SE; T-Stat; p-val
	#Interdep.
	y<-subset(d3, treat_nia==1 & is.na(support_123)==F)$support_123
		length(y)	#N
		sum(y==1)/(sum(y==1)+sum(y==0))	#Proportion Support
		t.test2( mean(y), mean(x), sd(y), sd(x), length(y), length(x))	#Diff; SE; T-Stat; p-val
	#Network
	y<-subset(d3, treat_nw==1 & is.na(support_123)==F)$support_123
		length(y)	#N
		sum(y==1)/(sum(y==1)+sum(y==0))	#Proportion Support
		t.test2( mean(y), mean(x), sd(y), sd(x), length(y), length(x))	#Diff; SE; T-Stat; p-val
	#Placebo
	y<-subset(d3, treat_plac==1 & is.na(support_123)==F)$support_123
		length(y)	#N
		sum(y==1)/(sum(y==1)+sum(y==0))	#Proportion Support
		t.test2( mean(y), mean(x), sd(y), sd(x), length(y), length(x))	#Diff; SE; T-Stat; p-val

	#DV: Support (Broad Defn) == 1 or 2 or 3 or 4
	#Null
	x<-subset(d3, treat_null==1 & is.na(support_1234)==F)$support_1234
		length(x)	#N
		sum(x==1)/(sum(x==1)+sum(x==0))	#Proportion Support
	#Direct
	y<-subset(d3, treat_oep==1 & is.na(support_1234)==F)$support_1234
		length(y)	#N
		sum(y==1)/(sum(y==1)+sum(y==0))	#Proportion Support
		t.test2( mean(y), mean(x), sd(y), sd(x), length(y), length(x))	#Diff; SE; T-Stat; p-val
	#Interdep.
	y<-subset(d3, treat_nia==1 & is.na(support_1234)==F)$support_1234
		length(y)	#N
		sum(y==1)/(sum(y==1)+sum(y==0))	#Proportion Support
		t.test2( mean(y), mean(x), sd(y), sd(x), length(y), length(x))	#Diff; SE; T-Stat; p-val
	#Network
	y<-subset(d3, treat_nw==1 & is.na(support_1234)==F)$support_1234
		length(y)	#N
		sum(y==1)/(sum(y==1)+sum(y==0))	#Proportion Support
		t.test2( mean(y), mean(x), sd(y), sd(x), length(y), length(x))	#Diff; SE; T-Stat; p-val
	#Placebo
	y<-subset(d3, treat_plac==1 & is.na(support_1234)==F)$support_1234
		length(y)	#N
		sum(y==1)/(sum(y==1)+sum(y==0))	#Proportion Support
		t.test2( mean(y), mean(x), sd(y), sd(x), length(y), length(x))	#Diff; SE; T-Stat; p-val

##############
#Table B.5 and
#Table B.6
##############

nrow(data)	#1049

	#Table B.5, Model (base)
		logit1 <-logit <- glm(y_12 ~ treat_oep + treat_nia + treat_nw + treat_plac,
		                      data = data, family = binomial)
		summary(logit)
		Direct<- data.frame()
		newdata = data.frame(treat_oep = 1,	treat_nia=0, treat_nw=0,treat_plac=0)
		Direct<- predict(logit, newdata, type="response")
		newdata = data.frame(treat_oep = 0,	treat_nia=1, treat_nw=0,treat_plac=0)
		Interdependence<- data.frame()
		Interdependence<- predict(logit, newdata, type="response")
		newdata = data.frame(treat_oep = 0,	treat_nia=0, treat_nw=1,treat_plac=0)
		Network<- data.frame()
		Network<- predict(logit, newdata, type="response")
		newdata = data.frame(treat_oep = 0,	treat_nia=0, treat_nw=0,treat_plac=1)
		Placebo<- data.frame()
		Placebo<- predict(logit, newdata, type="response")
		base<- data.frame()
		base<- rbind(Direct, Interdependence,Network,Placebo)
		
	#Table B.5, Model (1)
	logit_b5_1 <-logit <- glm(y_12 ~ treat_oep + treat_nia + treat_nw + treat_plac
			+ female 
			+ female*treat_oep
			+ female*treat_nia
			+ female*treat_nw
			+ female*treat_plac,
			data = data, family = binomial)
		summary(logit)
		#Table B5, Column
	femD<- data.frame()
	newdata = data.frame(female = 1, treat_oep = 1,	treat_nia=0, treat_nw=0,treat_plac=0)
	femD<- predict(logit, newdata, type="response")
	newdata = data.frame(female = 1, treat_oep = 0,	treat_nia=1, treat_nw=0,treat_plac=0)
	femI<- data.frame()
	femI<- predict(logit, newdata, type="response")
	newdata = data.frame(female = 1, treat_oep = 0,	treat_nia=0, treat_nw=1,treat_plac=0)
	femN<- data.frame()
	femN<- predict(logit, newdata, type="response")
	newdata = data.frame(female = 1, treat_oep = 0,	treat_nia=0, treat_nw=0,treat_plac=1)
	femP<- data.frame()
	femP<- predict(logit, newdata, type="response")
	fem<- data.frame()
	fem<- rbind(femD, femI,femN,femP)
	
	nonfemD<- data.frame()
	newdata = data.frame(female = 0, treat_oep = 1,	treat_nia=0, treat_nw=0,treat_plac=0)
	nonfemD<- predict(logit, newdata, type="response")
	newdata = data.frame(female = 0, treat_oep = 0,	treat_nia=1, treat_nw=0,treat_plac=0)
	nonfemI<- data.frame()
	nonfemI<- predict(logit, newdata, type="response")
	newdata = data.frame(female = 0, treat_oep = 0,	treat_nia=0, treat_nw=1,treat_plac=0)
	nonfemN<- data.frame()
	nonfemN<- predict(logit, newdata, type="response")
	newdata = data.frame(female = 0, treat_oep = 0,	treat_nia=0, treat_nw=0,treat_plac=1)
	nonfemP<- data.frame()
	nonfemP<- predict(logit, newdata, type="response")
	nonfem<- data.frame()
	nonfem<- rbind(nonfemD, nonfemI, nonfemN, nonfemP)
	
	
	#Table B.5, Model (2)
		(medage <- as.numeric(summary(data$age)[3]))
	data$agelow[data$age<medage] <-1
	data$agelow[data$age>=medage] <-0
		sum(data$agelow==1, na.rm=T)
		sum(data$agelow==0, na.rm=T)

	logit_b5_2 <-logit <- glm(y_12 ~ treat_oep + treat_nia + treat_nw + treat_plac
			+ agelow 
			+ agelow*treat_oep
			+ agelow*treat_nia
			+ agelow*treat_nw
			+ agelow*treat_plac,
			data = data, family = binomial)
		summary(logit)

	agelowD<- data.frame()
	newdata = data.frame(agelow = 1, treat_oep = 1,	treat_nia=0, treat_nw=0,treat_plac=0)
	agelowD<- predict(logit, newdata, type="response")
	newdata = data.frame(agelow = 1, treat_oep = 0,	treat_nia=1, treat_nw=0,treat_plac=0)
	agelowI<- data.frame()
	agelowI<- predict(logit, newdata, type="response")
	newdata = data.frame(agelow = 1, treat_oep = 0,	treat_nia=0, treat_nw=1,treat_plac=0)
	agelowN<- data.frame()
	agelowN<- predict(logit, newdata, type="response")
	newdata = data.frame(agelow = 1, treat_oep = 0,	treat_nia=0, treat_nw=0,treat_plac=1)
	agelowP<- data.frame()
	agelowP<- predict(logit, newdata, type="response")
	agelow<- data.frame()
	agelow<- rbind(agelowD, agelowI,agelowN,agelowP)

	agehighD<- data.frame()
	newdata = data.frame(agelow = 0, treat_oep = 1,	treat_nia=0, treat_nw=0,treat_plac=0)
	agehighD<- predict(logit, newdata, type="response")
	newdata = data.frame(agelow = 0, treat_oep = 0,	treat_nia=1, treat_nw=0,treat_plac=0)
	agehighI<- data.frame()
	agehighI<- predict(logit, newdata, type="response")
	newdata = data.frame(agelow = 0, treat_oep = 0,	treat_nia=0, treat_nw=1,treat_plac=0)
	agehighN<- data.frame()
	agehighN<- predict(logit, newdata, type="response")
	newdata = data.frame(agelow = 0, treat_oep = 0,	treat_nia=0, treat_nw=0,treat_plac=1)
	agehighP<- data.frame()
	agehighP<- predict(logit, newdata, type="response")
	agehigh<- data.frame()
	agehigh<- rbind(agehighD, agehighI,agehighN,agehighP)

	
	
	#Table B.5, Model (3)
	logit_b5_3 <-logit <- glm(y_12 ~ treat_oep + treat_nia + treat_nw + treat_plac
			+ white 
			+ white*treat_oep
			+ white*treat_nia
			+ white*treat_nw
			+ white*treat_plac,
			data = data, family = binomial)
		summary(logit)
	
	whiteD<- data.frame()
	newdata = data.frame(white = 1, treat_oep = 1,	treat_nia=0, treat_nw=0,treat_plac=0)
	whiteD<- predict(logit, newdata, type="response")
	newdata = data.frame(white = 1, treat_oep = 0,	treat_nia=1, treat_nw=0,treat_plac=0)
	whiteI<- data.frame()
	whiteI<- predict(logit, newdata, type="response")
	newdata = data.frame(white = 1, treat_oep = 0,	treat_nia=0, treat_nw=1,treat_plac=0)
	whiteN<- data.frame()
	whiteN<- predict(logit, newdata, type="response")
	newdata = data.frame(white = 1, treat_oep = 0,	treat_nia=0, treat_nw=0,treat_plac=1)
	whiteP<- data.frame()
	whiteP<- predict(logit, newdata, type="response")
	white<- data.frame()
	white<- rbind(whiteD, whiteI,whiteN,whiteP)

	nonwhiteD<- data.frame()
	newdata = data.frame(white = 0, treat_oep = 1,	treat_nia=0, treat_nw=0,treat_plac=0)
	nonwhiteD<- predict(logit, newdata, type="response")
	newdata = data.frame(white = 0, treat_oep = 0,	treat_nia=1, treat_nw=0,treat_plac=0)
	nonwhiteI<- data.frame()
	nonwhiteI<- predict(logit, newdata, type="response")
	newdata = data.frame(white = 0, treat_oep = 0,	treat_nia=0, treat_nw=1,treat_plac=0)
	nonwhiteN<- data.frame()
	nonwhiteN<- predict(logit, newdata, type="response")
	newdata = data.frame(white = 0, treat_oep = 0,	treat_nia=0, treat_nw=0,treat_plac=1)
	nonwhiteP<- data.frame()
	nonwhiteP<- predict(logit, newdata, type="response")
	nonwhite<- data.frame()
	nonwhite<- rbind(nonwhiteD, nonwhiteI,nonwhiteN,nonwhiteP)
	
	
	#Table B.5, Model (4)
	logit_b5_4 <-logit <- glm(y_12 ~ treat_oep + treat_nia + treat_nw + treat_plac
			+ educ_ba 
			+ educ_ba*treat_oep
			+ educ_ba*treat_nia
			+ educ_ba*treat_nw
			+ educ_ba*treat_plac,
			data = data, family = binomial)
		summary(logit)

	#################
	# Table B.5 output
	#################
	library(stargazer)
	stargazer(logit_b5_1,logit_b5_2,logit_b5_3,logit_b5_4, style="qje",
			align=T, ci=F,  	dep.var.labels.include = FALSE,
			no.space=TRUE)


	baD<- data.frame()
	newdata = data.frame(educ_ba = 1, treat_oep = 1,	treat_nia=0, treat_nw=0,treat_plac=0)
	baD<- predict(logit, newdata, type="response")
	newdata = data.frame(educ_ba = 1, treat_oep = 0,	treat_nia=1, treat_nw=0,treat_plac=0)
	baI<- data.frame()
	baI<- predict(logit, newdata, type="response")
	newdata = data.frame(educ_ba = 1, treat_oep = 0,	treat_nia=0, treat_nw=1,treat_plac=0)
	baN<- data.frame()
	baN<- predict(logit, newdata, type="response")
	newdata = data.frame(educ_ba = 1, treat_oep = 0,	treat_nia=0, treat_nw=0,treat_plac=1)
	baP<- data.frame()
	baP<- predict(logit, newdata, type="response")
	ba<- data.frame()
	ba<- rbind(baD, baI,baN,baP)
	
	nobaD<- data.frame()
	newdata = data.frame(educ_ba = 0, treat_oep = 1,	treat_nia=0, treat_nw=0,treat_plac=0)
	nobaD<- predict(logit, newdata, type="response")
	newdata = data.frame(educ_ba = 0, treat_oep = 0,	treat_nia=1, treat_nw=0,treat_plac=0)
	nobaI<- data.frame()
	nobaI<- predict(logit, newdata, type="response")
	newdata = data.frame(educ_ba = 0, treat_oep = 0,	treat_nia=0, treat_nw=1,treat_plac=0)
	nobaN<- data.frame()
	nobaN<- predict(logit, newdata, type="response")
	newdata = data.frame(educ_ba = 0, treat_oep = 0,	treat_nia=0, treat_nw=0,treat_plac=1)
	nobaP<- data.frame()
	nobaP<- predict(logit, newdata, type="response")
	noba<- data.frame()
	noba<- rbind(nobaD, nobaI,nobaN,nobaP)

	#################
	# Table B.6 output
	#################

	tableb6<- cbind(base,fem,nonfem,agelow,agehigh,white,nonwhite,ba,noba)	
	colnames(tableb6)<- c("Baseline", "Fem.","Non-fem.","Age low", "Age high", "White", "Non-white", "Yes BA", "No BA")
	
	install.packages("xtable")
	library(xtable)
	xtable(tableb6, digits=3)


##################
### FIGURE B.1
##################

	##################
	###"The Respondents"
	#	Re: Survey Respondents
	##################
	#"We recruited 1293 survey respondents..."
	nrow(d)	#1293

	#"1159 respondents completed the survey in a median time of 16 minutes"
			### Identify subset of respondents that (i) complete the survey and (ii) that provided a DV response
			#v10 = 'finished' == completed survey
			#consent_tx 	== gave consent
			#support_12		== provided value for DV
	nrow(d2<-subset(d, v10==1 & consent_tx=="Yes" & is.na(support_12)==F))
					#n = 1159 completion
	summary(d2$minutes)	#Median completion time 16 minutes

	#We limited our sample to respondents [between 8 and 32 minutes]
	#which exclud[es] the fastest 1.8% and slowest 7.7% of survey takers
		#Left-hand cutpoint at no less than 1/2 median response time:
		x <-8			#8 minutes
		nrow(subset(d2,minutes<x)) 		#drop 21 observations			
		nrow(subset(d2,minutes<x))/nrow(d2) #== 1.81% completed observations		

		# Right hand cutpoint at 32 minutes (2x median value)
		x <- 32			#32 minutes
		nrow(subset(d2,minutes>x)) 		#drop 89 observations			
		nrow(subset(d2,minutes>x))/nrow(d2) #== 7.68% completed observations

	### Create barplot to justify this selection
			#PREP THE BARPLOTS
				(vals <- sort(unique(d2$minutes))); summary(vals);length(vals) #78
				x <- 0:415;	count <- rep(0, length(x)); h <- as.data.frame(cbind(x,count))
					head(h)
				for(i in 1:length(vals)){h$count[h$x==vals[i]] <- sum(d2$minutes==vals[i])}
				B <- h$count	
			# SHOW DATA FOR 1-HOUR (60 mins)
			E <- B[1:61] 	#take values from t==0 to t==60 only
			sum(E)			#1126
				sum(B[62:411])	#32
	#####
	#Fig B.1
	#####
		#pdf("CompletionTimeDist60mins.pdf")
			(mids <- barplot(E, xlab="",ylab="",axes=F))
			axis(1, at=c(mids[1],mids[11],mids[21],mids[31],mids[41],mids[51],mids[61]),
					labels=c("0","10","20","30","40","50","60"))
			axis(2, las=2)
			abline(v=c((mids[8]+.5),mids[32+1]+.5), lty=2)
			mtext("Respondent Completion Time Distribution", side=3,cex=1.15, font=2, line=1.5)
				#mtext("Respondent Completion Time Distribution\n 0 to 60 minutes", side=3,cex=1.15, font=2, line=2)
				#mtext("(n=1126 displayed of 1159 observations)", side=3,cex=1.00, font=1, line=1.5)
				#mtext("33 obs > 60 minutes)", side=3,cex=1.00, font=1, line=0.5)
			mtext("minutes to completion", side=1,cex=1.15, font=1, line=2.5)
			mtext("number of respondents", side=2,cex=1.15, font=1, line=3)
		#dev.off()

######
#Figure B.2 (left)
######
		#setwd("C:/Users/ffung/Dropbox/Financial Regulations_ChaudoinWilf/Drafts_CW")
		#pdf("ModeratorDistribution_FR.pdf")
			(mids <- barplot(table(d3$frsum_red), xlab="",ylab="",axes=F))
			axis(1,at=c(0,mids[1],mids[2],mids[3],mids[4],5),c("","","","","",""))
			axis(2, las=2)
				summary(d3$frsum_red)
			mtext("Moderator Distribution\n Folk Realism", side=3,cex=1.15, font=2, line=1.5)
				#mtext("Respondent Completion Time Distribution\n 0 to 60 minutes", side=3,cex=1.15, font=2, line=2)
				#mtext("(n=1126 displayed of 1159 observations)", side=3,cex=1.00, font=1, line=1.5)
				#mtext("33 obs > 60 minutes)", side=3,cex=1.00, font=1, line=0.5)
			mtext("Folk Realism Response Count", side=1,cex=1.15, font=1, line=2.5)
			mtext("number of respondents", side=2,cex=1.15, font=1, line=3)
			abline(v=mids[1]+.6,lty=2)
			abline(v=mids[1]+.6,lty=2)
		#dev.off()

######
#Figure B.2 (right)
######
		#pdf("ModeratorDistribution_Eth.pdf")
		(mids<-barplot(table(d3$ethno_diff),xlab="", names.arg=F, axes=F))
			d3$ethno
			cbind(sort(unique(d3$ethno_diff)),1:nrow(mids))	#median value (-1.22) == mids[74])
				nrow(mids)			#171
				min(mids);max(mids)	#min = 0.7; max = 204.7
			axis(1,at=c(-10,mids[1],(mids[74]+.5), mids[172],250),c("","-6.25","-1.22","3.36","" ))#,las=2)
			axis(2,las=2)
			abline(v=(mids[74]+.5),lty=2)
				mtext("Ethnocentrism Measure", side=1,cex=1.15, font=1, line=2.5)
			mtext("number of respondents", side=2,cex=1.15, font=1, line=3)
			mtext("Moderator Distribution\n Ethnocentrism", side=3,cex=1.15, font=2, line=1.5)
		#dev.off()

###
# Figure B.3(a) (left)
###
	d<-d3
	attach(d)
	# Summarizes the number 1of respondents that received each treatment
	nbyti<- as.matrix(aggregate(support_12 ~ treatment , data=d, FUN=function(x) sum(!is.na(x)) ))
		n1<-nbyti[1,2]
		n2<-nbyti[2,2]
		n3<-nbyti[3,2]
		n4<-nbyti[4,2]
		n5<-nbyti[5,2]
	# Repeats that step, for the number of respondents by treatment who approved
	d2 <- d[ which(d$support_12 =="1"), ]
	abyti<- as.matrix(aggregate(support_12 ~ treatment, data=d2, FUN=function(x) sum(!is.na(x)) ))
		a1<-abyti[1,2]
		a2<-abyti[2,2]
		a3<-abyti[3,2]
		a4<-abyti[4,2]
		a5<-abyti[5,2]
	# Number of draws from the Beta posterior
		n<-5000
	# Draws from the posterior for each treatment
		pi_1 <- as.matrix(rbeta(n, a1+0.5, n1-a1+0.5))
		pi_2 <- as.matrix(rbeta(n, a2+0.5, n2-a2+0.5))
		pi_3 <- as.matrix(rbeta(n, a3+0.5, n3-a3+0.5))
		pi_4 <- as.matrix(rbeta(n, a4+0.5, n4-a4+0.5))
		pi_5 <- as.matrix(rbeta(n, a5+0.5, n5-a5+0.5))
	# Treatment indicators to make the data for the plot cleaner
			#TREATMENT KEY / CODING
				#	1 = OEP
				#	2 = nia / interdependence
				#	3 = NW / network
				#	4 = null
				#	5 = placebo
		t1 <- as.matrix(rep(1,n, nrow=n, ncol=1))
		t2 <- as.matrix(rep(2,n, nrow=n, ncol=1))
		t3 <- as.matrix(rep(3,n, nrow=n, ncol=1))
		t4 <- as.matrix(rep(4,n, nrow=n, ncol=1))
		t5 <- as.matrix(rep(5,n, nrow=n, ncol=1))
	tis <- as.matrix(rbind(t1,t2,t3,t4,t5))
	pis <- as.matrix(rbind(pi_1,pi_2,pi_3,pi_4,pi_5))

	betas <- as.data.frame(cbind(tis,pis))
	betas$treatmentname = factor(betas$V1, labels = c("Direct","Interdep.","Network","Null","Plac."))

	
	#pdf("PercSupport_12_ByTreatment_BetaEsts_Appendix.pdf")
	lineplot.CI(x.factor = betas$treatmentname, response = betas$V2, type="p", ylab="",xlab="",axes=F, font=2,
			data = betas, main = "", #xlab = "Treatment Group", ylab = "% Supporting Regulation (12)",
			ylim=c(0,1),
			ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)) )
		mtext("Percent Supporting Regulation \n Narrow Definition of Support", side=3,cex=1.15, font=2, line=1.5)
		axis(2, las=2)
		mtext("Treatment Group", side=1,cex=1.15, font=1, line=2.5)
	#dev.off()

###
# Figure B.3(b) (middle)
###
	###
	# Support 123, Betas
	###
	# Summarizes the number of respondents that received each treatment
	nbyti<- as.matrix(aggregate(support_123 ~ treatment , data=d, FUN=function(x) sum(!is.na(x)) ))
		n1<-nbyti[1,2]
		n2<-nbyti[2,2]
		n3<-nbyti[3,2]
		n4<-nbyti[4,2]
		n5<-nbyti[5,2]
	# Repeats that step, for the number of respondents by treatment who approved
	d2 <- d[ which(d$support_123 =="1"), ]
	abyti<- as.matrix(aggregate(support_123 ~ treatment, data=d2, FUN=function(x) sum(!is.na(x)) ))
		a1<-abyti[1,2]
		a2<-abyti[2,2]
		a3<-abyti[3,2]
		a4<-abyti[4,2]
		a5<-abyti[5,2]
	# Number of draws from the Beta posterior
		n<-5000
	# Draws from the posterior for each treatment
		pi_1 <- as.matrix(rbeta(n, a1+0.5, n1-a1+0.5))
		pi_2 <- as.matrix(rbeta(n, a2+0.5, n2-a2+0.5))
		pi_3 <- as.matrix(rbeta(n, a3+0.5, n3-a3+0.5))
		pi_4 <- as.matrix(rbeta(n, a4+0.5, n4-a4+0.5))
		pi_5 <- as.matrix(rbeta(n, a5+0.5, n5-a5+0.5))
	# Treatment indicators to make the data for the plot cleaner
		t1 <- as.matrix(rep(1,n, nrow=n, ncol=1))
		t2 <- as.matrix(rep(2,n, nrow=n, ncol=1))
		t3 <- as.matrix(rep(3,n, nrow=n, ncol=1))
		t4 <- as.matrix(rep(4,n, nrow=n, ncol=1))
		t5 <- as.matrix(rep(5,n, nrow=n, ncol=1))
	tis <- as.matrix(rbind(t1,t2,t3,t4,t5))
	pis <- as.matrix(rbind(pi_1,pi_2,pi_3,pi_4,pi_5))

	betas <- as.data.frame(cbind(tis,pis))
	betas$treatmentname = factor(betas$V1, labels = c("Direct","Interdep.","Network","Null","Plac."))

	# Approval ratings by treatment group
	#pdf("PercSupport_123_ByTreatment_BetaEsts_Appendix.pdf")
	lineplot.CI(x.factor = betas$treatmentname, response = betas$V2, type="p", ylab="",xlab="",axes=F, font=2,
			data = betas, main = "", 
			ylim=c(0,1),
			ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)) )
		mtext("Percent Supporting Regulation \n Moderate Definition of Support", side=3,cex=1.15, font=2, line=1.5)
		axis(2, las=2)
		mtext("Treatment Group", side=1,cex=1.15, font=1, line=2.5)
	#dev.off()

###
# Figure B.3(c) (right)
###
	###
	# Support 1234
	###
	# Summarizes the number of respondents that received each treatment
	nbyti<- as.matrix(aggregate(support_1234 ~ treatment , data=d, FUN=function(x) sum(!is.na(x)) ))
		n1<-nbyti[1,2]
		n2<-nbyti[2,2]
		n3<-nbyti[3,2]
		n4<-nbyti[4,2]
		n5<-nbyti[5,2]
	# Repeats that step, for the number of respondents by treatment who approved
	d2 <- d[ which(d$support_1234 =="1"), ]
	abyti<- as.matrix(aggregate(support_1234 ~ treatment, data=d2, FUN=function(x) sum(!is.na(x)) ))
		a1<-abyti[1,2]
		a2<-abyti[2,2]
		a3<-abyti[3,2]
		a4<-abyti[4,2]
		a5<-abyti[5,2]
	# Number of draws from the Beta posterior
		n<-5000
	# Draws from the posterior for each treatment
	pi_1 <- as.matrix(rbeta(n, a1+0.5, n1-a1+0.5))
	pi_2 <- as.matrix(rbeta(n, a2+0.5, n2-a2+0.5))
	pi_3 <- as.matrix(rbeta(n, a3+0.5, n3-a3+0.5))
	pi_4 <- as.matrix(rbeta(n, a4+0.5, n4-a4+0.5))
	pi_5 <- as.matrix(rbeta(n, a5+0.5, n5-a5+0.5))

	# Treatment indicators to make the data for the plot cleaner
	t1 <- as.matrix(rep(1,n, nrow=n, ncol=1))
	t2 <- as.matrix(rep(2,n, nrow=n, ncol=1))
	t3 <- as.matrix(rep(3,n, nrow=n, ncol=1))
	t4 <- as.matrix(rep(4,n, nrow=n, ncol=1))
	t5 <- as.matrix(rep(5,n, nrow=n, ncol=1))

	tis <- as.matrix(rbind(t1,t2,t3,t4,t5))
	pis <- as.matrix(rbind(pi_1,pi_2,pi_3,pi_4,pi_5))

	betas <- as.data.frame(cbind(tis,pis))
	betas$treatmentname = factor(betas$V1, labels = c("Direct","Interdep.","Network","Null","Plac."))
	###
	# Figure B.3(c) (right)
	###
	# Approval ratings by treatment group
	#pdf("PercSupport_1234_ByTreatment_BetaEsts_Appendix.pdf")
	lineplot.CI(x.factor = betas$treatmentname, response = betas$V2, type="p", ylab="",xlab="",axes=F, font=2,
			data = betas, main = "", 
			ylim=c(0,1),
			ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)) )
		mtext("Percent Supporting Regulation \n Broad Definition of Support", side=3,cex=1.15, font=2, line=1.5)
		axis(2, las=2)
		mtext("Treatment Group", side=1,cex=1.15, font=1, line=2.5)
	#dev.off()


##############
#Figure B.4 
# High and low values of Folk Realism, Ethnocentrism
#############

	d3$ansamp_ethnoh_ind <- ifelse(is.na(d3$ethno_diff)==F,ifelse(d3$ethno_diff>summary(d3$ethno_diff)[3],1,0),NA)
	d <- d3
		d_frh <- d[ which(d$frsum_red != 0), ]
			nrow(d_frh)	#677
		d_frl <- d[ which(d$frsum_red == 0), ]
			nrow(d_frl)	#363
		d_ethh <- d[ which(d$ansamp_ethnoh_ind == 1), ]
			nrow(d_ethh)	#448
		d_ethl <- d[ which(d$ansamp_ethnoh_ind == 0), ]
			nrow(d_ethl)	#465
	###
	# Support 12
	# Treatment pictures for H/L, H/L
	###

	# Summarizes the number of respondents that received each treatment, by h/l
	nbyti_frh <- as.matrix(aggregate(support_12 ~ treatment , data=d_frh, FUN=function(x) sum(!is.na(x)) ))
	nbyti_frl <- as.matrix(aggregate(support_12 ~ treatment , data=d_frl, FUN=function(x) sum(!is.na(x)) ))
	nbyti_ethh <- as.matrix(aggregate(support_12 ~ treatment , data=d_ethh, FUN=function(x) sum(!is.na(x)) ))
	nbyti_ethl <- as.matrix(aggregate(support_12 ~ treatment , data=d_ethl, FUN=function(x) sum(!is.na(x)) ))

	n1_frh<-nbyti_frh[1,2]
	n2_frh<-nbyti_frh[2,2]
	n3_frh<-nbyti_frh[3,2]
	n4_frh<-nbyti_frh[4,2]
	n5_frh<-nbyti_frh[5,2]

	n1_frl<-nbyti_frl[1,2]
	n2_frl<-nbyti_frl[2,2]
	n3_frl<-nbyti_frl[3,2]
	n4_frl<-nbyti_frl[4,2]
	n5_frl<-nbyti_frl[5,2]

	n1_ethh<-nbyti_ethh[1,2]
	n2_ethh<-nbyti_ethh[2,2]
	n3_ethh<-nbyti_ethh[3,2]
	n4_ethh<-nbyti_ethh[4,2]
	n5_ethh<-nbyti_ethh[5,2]

	n1_ethl<-nbyti_ethl[1,2]
	n2_ethl<-nbyti_ethl[2,2]
	n3_ethl<-nbyti_ethl[3,2]
	n4_ethl<-nbyti_ethl[4,2]
	n5_ethl<-nbyti_ethl[5,2]


	# Repeats that step, for the number of respondents by treatment who approved
	d2_frh <- d_frh[ which(d_frh$support_12 =="1"), ]
	d2_frl <- d_frl[ which(d_frl$support_12 =="1"), ]
	d2_ethh <- d_ethh[ which(d_ethh$support_12 =="1"), ]
	d2_ethl <- d_ethl[ which(d_ethl$support_12 =="1"), ]

	abyti_frh<- as.matrix(aggregate(support_12 ~ treatment, data=d2_frh, FUN=function(x) sum(!is.na(x)) ))
	abyti_frl<- as.matrix(aggregate(support_12 ~ treatment, data=d2_frl, FUN=function(x) sum(!is.na(x)) ))
	abyti_ethh<- as.matrix(aggregate(support_12 ~ treatment, data=d2_ethh, FUN=function(x) sum(!is.na(x)) ))
	abyti_ethl<- as.matrix(aggregate(support_12 ~ treatment, data=d2_ethl, FUN=function(x) sum(!is.na(x)) ))

	a1_frh<-abyti_frh[1,2]
	a2_frh<-abyti_frh[2,2]
	a3_frh<-abyti_frh[3,2]
	a4_frh<-abyti_frh[4,2]
	a5_frh<-abyti_frh[5,2]

	a1_frl<-abyti_frl[1,2]
	a2_frl<-abyti_frl[2,2]
	a3_frl<-abyti_frl[3,2]
	a4_frl<-abyti_frl[4,2]
	a5_frl<-abyti_frl[5,2]

	a1_ethh<-abyti_ethh[1,2]
	a2_ethh<-abyti_ethh[2,2]
	a3_ethh<-abyti_ethh[3,2]
	a4_ethh<-abyti_ethh[4,2]
	a5_ethh<-abyti_ethh[5,2]

	a1_ethl<-abyti_ethl[1,2]
	a2_ethl<-abyti_ethl[2,2]
	a3_ethl<-abyti_ethl[3,2]
	a4_ethl<-abyti_ethl[4,2]
	a5_ethl<-abyti_ethl[5,2]
	
	# Number of draws from the Beta posterior
	n<-5000

	# Draws from the posterior for each treatment
	pi_1_frh <- as.matrix(rbeta(n, a1_frh+0.5, n1_frh-a1_frh+0.5))
	pi_2_frh <- as.matrix(rbeta(n, a2_frh+0.5, n2_frh-a2_frh+0.5))
	pi_3_frh <- as.matrix(rbeta(n, a3_frh+0.5, n3_frh-a3_frh+0.5))
	pi_4_frh <- as.matrix(rbeta(n, a4_frh+0.5, n4_frh-a4_frh+0.5))
	pi_5_frh <- as.matrix(rbeta(n, a5_frh+0.5, n5_frh-a5_frh+0.5))

	pi_1_frl <- as.matrix(rbeta(n, a1_frl+0.5, n1_frl-a1_frl+0.5))
	pi_2_frl <- as.matrix(rbeta(n, a2_frl+0.5, n2_frl-a2_frl+0.5))
	pi_3_frl <- as.matrix(rbeta(n, a3_frl+0.5, n3_frl-a3_frl+0.5))
	pi_4_frl <- as.matrix(rbeta(n, a4_frl+0.5, n4_frl-a4_frl+0.5))
	pi_5_frl <- as.matrix(rbeta(n, a5_frl+0.5, n5_frl-a5_frl+0.5))

	pi_1_ethh <- as.matrix(rbeta(n, a1_ethh+0.5, n1_ethh-a1_ethh+0.5))
	pi_2_ethh <- as.matrix(rbeta(n, a2_ethh+0.5, n2_ethh-a2_ethh+0.5))
	pi_3_ethh <- as.matrix(rbeta(n, a3_ethh+0.5, n3_ethh-a3_ethh+0.5))
	pi_4_ethh <- as.matrix(rbeta(n, a4_ethh+0.5, n4_ethh-a4_ethh+0.5))
	pi_5_ethh <- as.matrix(rbeta(n, a5_ethh+0.5, n5_ethh-a5_ethh+0.5))

	pi_1_ethl <- as.matrix(rbeta(n, a1_ethl+0.5, n1_ethl-a1_ethl+0.5))
	pi_2_ethl <- as.matrix(rbeta(n, a2_ethl+0.5, n2_ethl-a2_ethl+0.5))
	pi_3_ethl <- as.matrix(rbeta(n, a3_ethl+0.5, n3_ethl-a3_ethl+0.5))
	pi_4_ethl <- as.matrix(rbeta(n, a4_ethl+0.5, n4_ethl-a4_ethl+0.5))
	pi_5_ethl <- as.matrix(rbeta(n, a5_ethl+0.5, n5_ethl-a5_ethl+0.5))


	# Treatment indicators to make the data for the plot cleaner
	t1 <- as.matrix(rep(1,n, nrow=n, ncol=1))
	t2 <- as.matrix(rep(2,n, nrow=n, ncol=1))
	t3 <- as.matrix(rep(3,n, nrow=n, ncol=1))
	t4 <- as.matrix(rep(4,n, nrow=n, ncol=1))
	t5 <- as.matrix(rep(5,n, nrow=n, ncol=1))

	tis <- as.matrix(rbind(t1,t2,t3,t4,t5))

	pis_frh <- as.matrix(rbind(pi_1_frh,pi_2_frh,pi_3_frh,pi_4_frh,pi_5_frh))
	pis_frl <- as.matrix(rbind(pi_1_frl,pi_2_frl,pi_3_frl,pi_4_frl,pi_5_frl))
	pis_ethh <- as.matrix(rbind(pi_1_ethh,pi_2_ethh,pi_3_ethh,pi_4_ethh,pi_5_ethh))
	pis_ethl <- as.matrix(rbind(pi_1_ethl,pi_2_ethl,pi_3_ethl,pi_4_ethl,pi_5_ethl))

	betas_frh <- as.data.frame(cbind(tis,pis_frh))
	betas_frl <- as.data.frame(cbind(tis,pis_frl))
	betas_ethh <- as.data.frame(cbind(tis,pis_ethh))
	betas_ethl <- as.data.frame(cbind(tis,pis_ethl))

	betas_frh$treatmentname = factor(betas_frh$V1, labels = c("Direct","Interdep.","Network","Null","Plac."))
	betas_frl$treatmentname = factor(betas_frl$V1, labels = c("Direct","Interdep.","Network","Null","Plac."))
	betas_ethh$treatmentname = factor(betas_ethh$V1, labels = c("Direct","Interdep.","Network","Null","Plac."))
	betas_ethl$treatmentname = factor(betas_ethl$V1, labels = c("Direct","Interdep.","Network","Null","Plac."))


	###
	# Figure B.4 (1 of 4 - top left)
	###
	#pdf("PercSupport_12_Mod_FolkRealismHigh_BetaEsts.pdf")
		# Approval ratings by treatment group, frh
		lineplot.CI(x.factor = betas_frh$treatmentname, response = betas_frh$V2, type="p",
	data = betas_frh, #main = "_frh", xlab = "Treatment Group", ylab = "% Supporting Regulation (Medium)",
	ylim=c(0,1),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)) )
	mtext("Folk Realism High", side=3,cex=2, font=2, line=1)
	#dev.off()

	###
	# Figure B.4 (2 of 4 - top right)
	###
	#pdf("PercSupport_12_Mod_FolkRealismLow_BetaEsts.pdf")
		# Approval ratings by treatment group, frl
		lineplot.CI(x.factor = betas_frl$treatmentname, response = betas_frl$V2, type="p",
	data = betas_frl, #main = "_frl", xlab = "Treatment Group", ylab = "% Supporting Regulation (Medium)",
	ylim=c(0,1),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)) )
	mtext("Folk Realism Low", side=3,cex=2, font=2, line=1)
	#dev.off()

	###
	# Figure B.4 (3 of 4 - bottom left)
	###
	#pdf("PercSupport_12_Mod_EthnoHigh_BetaEsts.pdf")
		# Approval ratings by treatment group, ethh
		lineplot.CI(x.factor = betas_ethh$treatmentname, response = betas_ethh$V2, type="p",
	data = betas_ethh, #main = "_ethh", xlab = "Treatment Group", ylab = "% Supporting Regulation (Medium)",
	ylim=c(0,1),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)) )
	mtext("Ethnocentrism High", side=3,cex=2, font=2, line=1)
	#dev.off()

	###
	# Figure B.4 (4 of 4 - bottom right)
	###
	#pdf("PercSupport_12_Mod_EthnoLow_BetaEsts.pdf")
		# Approval ratings by treatment group, ethl
	lineplot.CI(x.factor = betas_ethl$treatmentname, response = betas_ethl$V2, type="p",
	data = betas_ethl, #main = "_ethl", xlab = "Treatment Group", ylab = "% Supporting Regulation (Medium)",
	ylim=c(0,1),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)) )
		mtext("Ethnocentrism Low", side=3,cex=2, font=2, line=1)
	#mtext("(n=)", side=3,cex=1.25, font=2, line=0)
		#dev.print(postscript, "/Users/robertchaudoin/Dropbox/Financial Regulations/Drafts_CW/support_by_tmt12_tls_ethl.eps", horizontal=FALSE)
	#dev.off()

####################
#Figure 3
####################

	####################
	#Prepare for Figure 3
	####################

	d <- read.csv("ChaudoinWilf_FPA_data.csv")
	nrow(d3<-subset(d, v10==1 & consent_tx=="Yes" & is.na(support_12)==F & minutes>=8 & minutes<=32)) #n = 1049
	d<-d3

	# This is a function that makes any of the beta figures.
	#	dsubset is the subset (or full) part of the data you're using to make the figure
	make_beta_fig <- function(dsubset)	{
	# Summarizes the number of respondents that received each treatment
	nbyti<- as.matrix(aggregate(support_12 ~ treatment , data= dsubset, FUN=function(x) sum(!is.na(x)) ))
		n1<-nbyti[1,2]
		n2<-nbyti[2,2]
		n3<-nbyti[3,2]
		n4<-nbyti[4,2]
		n5<-nbyti[5,2]
	# Repeats that step, for the number of respondents by treatment who approved
	d2 <- dsubset[ which(dsubset $support_12 =="1"), ]
	abyti<- as.matrix(aggregate(support_12 ~ treatment, data=d2, FUN=function(x) sum(!is.na(x)) ))
		a1<-abyti[1,2]
		a2<-abyti[2,2]
		a3<-abyti[3,2]
		a4<-abyti[4,2]
		a5<-abyti[5,2]
	# Number of draws from the Beta posterior
		n<-5000
	# Draws from the posterior for each treatment
		pi_1 <- as.matrix(rbeta(n, a1+0.5, n1-a1+0.5))
		pi_2 <- as.matrix(rbeta(n, a2+0.5, n2-a2+0.5))
		pi_3 <- as.matrix(rbeta(n, a3+0.5, n3-a3+0.5))
		pi_4 <- as.matrix(rbeta(n, a4+0.5, n4-a4+0.5))
		pi_5 <- as.matrix(rbeta(n, a5+0.5, n5-a5+0.5))
	# Treatment indicators to make the data for the plot cleaner
		t1 <- as.matrix(rep(1,n, nrow=n, ncol=1))
		t2 <- as.matrix(rep(2,n, nrow=n, ncol=1))
		t3 <- as.matrix(rep(3,n, nrow=n, ncol=1))
		t4 <- as.matrix(rep(4,n, nrow=n, ncol=1))
		t5 <- as.matrix(rep(5,n, nrow=n, ncol=1))
	tis <- as.matrix(rbind(t1,t2,t3,t4,t5))
	pis <- as.matrix(rbind(pi_1,pi_2,pi_3,pi_4,pi_5))
	betas <- as.data.frame(cbind(tis,pis))
	betas$treatmentname = factor(betas$V1, labels = c("Direct","Interdep.","Network","Null","Plac."))
	print(nrow(betas))
	print(summary(betas$V2))
	print(names(betas))
	return(betas)
	}
# Making Beta figures with each subset of the data
##########
#Figure 3 (1 of 4 - top left)
##########
# FR High
plot.new()
frh <- d[ which(d$fr_123_ind == 1), ]
		#frh <- d[ which(d$fr_23_ind == 1), ]
betas_frh <- make_beta_fig(frh)
	#pdf("fromfunc_betas_frh_nia_null.pdf")
			#	pdf("/Users/robertchaudoin/Dropbox/Financial Regulations/Drafts_CW/fromfunc_betas_frh_nia_null.pdf")
betas_frh_nia_null <- betas_frh[ which(betas_frh$treatmentname == "Interdep." | betas_frh$treatmentname == "Null"), ]
		lineplot.CI(x.factor = betas_frh_nia_null$treatmentname, response = betas_frh_nia_null$V2, type="p",
			data = betas_frh_nia_null, xlab="",ylab="",#main = "Folk Realism: High", xlab = "Treatment Group", ylab = "% Supporting Regulation (12)",
			ylim=c(0.3,0.8),axes=F,
			ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)) )
	mtext("Folk Realists", side=3,cex=2, font=2, line=1)
	mtext("(High levels of folk realism)", side=3,cex=1.25, font=2, line=0)
		axis(2, las=2)
		mtext("Treatment Group", side=1,cex=1.15, font=1, line=2.5)
	#dev.off()
##########
#Figure 3 (2 of 4 - top right)
##########
# FR Low
frl <- d[ which(d$fr_123_ind == 0), ]
		#frl <- d[ which(d$fr_23_ind == 0), ]
betas_frl <- make_beta_fig(frl)
	#pdf("fromfunc_betas_frl_nia_null.pdf")
			#	pdf("/Users/robertchaudoin/Dropbox/Financial Regulations/Drafts_CW/fromfunc_betas_frl_nia_null.pdf")
betas_frl_nia_null <- betas_frl[ which(betas_frl$treatmentname == "Interdep." | betas_frl$treatmentname == "Null"), ]
		lineplot.CI(x.factor = betas_frl_nia_null$treatmentname, response = betas_frl_nia_null $V2, type="p",
			data = betas_frl_nia_null, xlab="",ylab="",#main = "Folk Realism: Low", xlab = "Treatment Group", ylab = "% Supporting Regulation (12)",
			ylim=c(0.3,0.8),axes=F,
			ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)) )
	mtext("Non-Folk Realists", side=3,cex=2, font=2, line=1)
	mtext("(Low levels of folk realism)", side=3,cex=1.25, font=2, line=0)
		axis(2, las=2)
		mtext("Treatment Group", side=1,cex=1.15, font=1, line=2.5)
	#dev.off()
##########
#Figure 3 (3 of 4 - bottom left)
##########
# Eth High
ethh <- d[ which(d$ansamp_ethnoh_ind == 1), ]
betas_ethh <- make_beta_fig(ethh)
	#pdf("fromfunc_betas_ethh_nia_null.pdf")
			#	pdf("/Users/robertchaudoin/Dropbox/Financial Regulations/Drafts_CW/fromfunc_betas_ethh_nia_null.pdf")
betas_ethh_nia_null <- betas_ethh[ which(betas_ethh$treatmentname == "Interdep." | betas_ethh$treatmentname == "Null"), ]
		lineplot.CI(x.factor = betas_ethh_nia_null$treatmentname, response = betas_ethh_nia_null$V2, type="p",
			data = betas_ethh_nia_null, xlab="",ylab="",# main = "Ethnocentrism: High", xlab = "Treatment Group", ylab = "% Supporting Regulation (12)",
			ylim=c(0.2,0.65),axes=F,
			ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)) )
	mtext("Ethnocentrist", side=3,cex=2, font=2, line=1)
	mtext("(High levels of ethnocentrism)", side=3,cex=1.25, font=2, line=0)
		axis(2, las=2)
		mtext("Treatment Group", side=1,cex=1.15, font=1, line=2.5)
	#dev.off()
##########
#Figure 3 (4 of 4 - bottom right)
##########
# Eth Low
ethl <- d[ which(d$ansamp_ethnoh_ind == 0), ]
betas_ethl <- make_beta_fig(ethl)
	#pdf("fromfunc_betas_ethl_nia_null.pdf")
			#	pdf("/Users/robertchaudoin/Dropbox/Financial Regulations/Drafts_CW/fromfunc_betas_ethl_nia_null.pdf")
betas_ethl_nia_null <- betas_ethl[ which(betas_ethl$treatmentname == "Interdep." | betas_ethl$treatmentname == "Null"), ]
		lineplot.CI(x.factor = betas_ethl_nia_null$treatmentname, response = betas_ethl_nia_null$V2, type="p",
			data = betas_ethl_nia_null, xlab="",ylab="",# main = "Ethnocentrism: Low", xlab = "Treatment Group", ylab = "% Supporting Regulation (12)",
			ylim=c(0.2,0.65),axes=F,
			ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)) )
	mtext("Non-Ethnocentrist", side=3,cex=2, font=2, line=1)
	mtext("(Low levels of ethnocentrism)", side=3,cex=1.25, font=2, line=0)
		axis(2, las=2)
		mtext("Treatment Group", side=1,cex=1.15, font=1, line=2.5)
	#dev.off()
